﻿#region << 版 本 注 释 >>
/*----------------------------------------------------------------
 * 版权所有 (c) 2022 北京超维景生物科技有限公司 保留所有权利。
 * 
 * 创建者：huangyang
 * 电子邮箱：huangyang@tvscope.cn
 * 创建时间：2023/3/9 14:56:23
 * 版本：V1.0.0
 * 描述：
 *
 * ----------------------------------------------------------------
 * 修改人：
 * 时间：
 * 修改说明：
 *
 * 版本：V1.0.1
 *----------------------------------------------------------------*/
#endregion << 版 本 注 释 >>

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace ImageK.Measure
{
    public class SplineFitter
    {
        private double[] y2;
        private static int EXTEND_BY = 7;
        private int extendBy;
        private float[] xpoints, ypoints;
        private int npoints;
        private int[] ixpoints, iypoints;

        public SplineFitter(int[] x, int[] y, int n)
        {
            // initSpline(x, y, n);
        }

        /** For closed curves: the first and last y value should be identical;
		 *	internally, a periodic continuation with a few will be used at both
		 *	ends */
        public SplineFitter(float[] x, float[] y, int n, bool closed)
        {
            // initSpline(x, y, n, closed);
        }

        public SplineFitter(float[] x, float[] y, int n)
        {
            // initSpline(x, y, n, false);
        }

		/** Given arrays of data points x[0..n-1] and y[0..n-1], computes the
		values of the second derivative at each of the data points
		y2[0..n-1] for use in the evalSpline() function. */
		private void initSpline(int[] x, int[] y, int n)
		{
			int i, k;
			double p, qn, sig, un;
			y2 = new double[n];  // cached
			double[] u = new double[n];
			for (i = 1; i < n - 1; i++)
			{
				// 888 chk for div by 0?
				sig = ((double)x[i] - x[i - 1]) / ((double)x[i + 1] - x[i - 1]);
				p = sig * y2[i - 1] + 2.0;
				y2[i] = (sig - 1.0) / p;
				u[i] = (((double)y[i + 1] - y[i]) / (x[i + 1] - x[i])) -
					   (((double)y[i] - y[i - 1]) / (x[i] - x[i - 1]));
				u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
			}
			qn = un = 0.0;
			y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
			for (k = n - 2; k >= 0; k--)
				y2[k] = y2[k] * y2[k + 1] + u[k];
			ixpoints = x;
			iypoints = y;
			npoints = n;
		}

		private void initSpline(float[] x, float[] y, int n, bool closed)
		{
			if (closed)
			{                   //add periodic continuation at both ends
				extendBy = EXTEND_BY;
				if (extendBy >= n)
					extendBy = n - 1;
				int n2 = n + 2 * extendBy;
				float[] xx = new float[n2];
				float[] yy = new float[n2];
				for (int i = 0; i < extendBy; i++)
				{
					xx[i] = x[n - (extendBy - i + 1)] - x[n - 1];
					yy[i] = y[n - (extendBy - i + 1)];
				}
				for (int i = extendBy; i < extendBy + n; i++)
				{
					xx[i] = x[i - extendBy];
					yy[i] = y[i - extendBy];
				}
				for (int i = extendBy + n; i < n2; i++)
				{
					xx[i] = x[i + 1 - (extendBy + n)] - x[0] + x[n - 1];
					yy[i] = y[i + 1 - (extendBy + n)];
				}
				n = n2;
				x = xx;
				y = yy;
			}

            {
                int i, k;
                double p, qn, sig, un;
                y2 = new double[n]; // cached
                double[] u = new double[n];
                for (i = 1; i < n - 1; i++)
                {
                    // 888 chk for div by 0?
                    sig = ((double)x[i] - x[i - 1]) / ((double)x[i + 1] - x[i - 1]);
                    p = sig * y2[i - 1] + 2.0;
                    y2[i] = (sig - 1.0) / p;
                    u[i] = (((double)y[i + 1] - y[i]) / (x[i + 1] - x[i])) -
                           (((double)y[i] - y[i - 1]) / (x[i] - x[i - 1]));
                    u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
                }

                qn = un = 0.0;
                y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
                for (k = n - 2; k >= 0; k--)
                    y2[k] = y2[k] * y2[k + 1] + u[k];
                xpoints = x;
                ypoints = y;
                npoints = n;
            }
        }

		/** Evalutes spline function at given point */
		public double evalSpline(double xp)
		{
			if (xpoints != null)
				return evalSpline(xpoints, ypoints, npoints, xp);
			else
				return evalSpline(ixpoints, iypoints, npoints, xp);
		}

		public double evalSpline(int[] x, int[] y, int n, double xp)
		{
			int klo, khi, k;
			double h, b, a;
			klo = 0;
			khi = n - 1;
			while (khi - klo > 1)
			{
				k = (khi + klo) >> 1;
				if (x[k] > xp) khi = k;
				else klo = k;
			}
			h = x[khi] - x[klo];
			/* orig code */
			/* if (h==0.0) FatalError("bad xvalues in splint\n"); */
			if (h == 0.0) return (0.0);  /* arbitr ret for now */
			a = (x[khi] - xp) / h;
			b = (xp - x[klo]) / h;
			// should have better err checking
			if (y2 == null) return (0.0);
			return (a * y[klo] + b * y[khi] + ((a * a * a - a) * y2[klo] + (b * b * b - b) * y2[khi]) * (h * h) / 6.0);
		}

		public double evalSpline(float[] x, float[] y, int n, double xp)
		{
			int klo, khi, k;
			double h, b, a;
			klo = 0;
			khi = n - 1;
			while (khi - klo > 1)
			{
				k = (khi + klo) >> 1;
				if (x[k] > xp)
					khi = k;
				else
					klo = k;
			}
			h = x[khi] - x[klo];
			/* orig code */
			/* if (h==0.0) FatalError("bad xvalues in splint\n"); */
			if (h == 0.0)
				return (0.0);  /* arbitr ret for now */
			a = (x[khi] - xp) / h;
			b = (xp - x[klo]) / h;
			// should have better err checking
			if (y2 == null)
				return (0.0);
			return (a * y[klo] + b * y[khi] + ((a * a * a - a) * y2[klo] + (b * b * b - b) * y2[khi]) * (h * h) / 6.0);
		}

	}
}
